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Abstract. The solutions of the Balitsky-Kovchegov evolution equations are studied numerically and com- 
pared with known analytical estimations. The rapidity and nuclear size dependences of the saturation scale 
are obtained for the cases of fixed and running coupling constant. These same dependences are studied in 
experimental data, on lepton-nucleus, deuteron-nucleus and nucleus-nucleus collisions, through geometric 
scaling and compared with the theoretical calculations. 
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1 Introduction 

The cleanest experimental information about parton dis- 
tributions comes from deep inelastic scattering experi- 
ments. At high energy, they can be described by the QCD 
dipole model [1] , which expresses the cross section of the 
virtual photon, emitted by the lepton, on the hadron h 
(proton or nucleus) as 



2 BK equation 

The BK equation [2] describes the rapidity Y = ln(s/so) = 
ln(a;o/x) evolution of the scattering probability N(x, y, Y) 
of a qq dipole with an hadronic target. When considering 
an homogeneous target with radius much larger than the 
size of any dipole, the dependence on impact parameter 
can be neglected and the equation reads 



dr 



dz\^ L (Q 2 ,v,z)\ 2 ai p (v,x). (1) 9N(r,Y) 



d 2 z 



K{r,r u r 2 ) N{r\,Y) 



Here, \Pt.l are the transverse and longitudinal wave func- 
tions for the splitting of 7* into a qq dipole of transverse 
size r with light-cone fractions z and (1 — z). The dipole 
cross section, cr^ ip (r,x), can be written as an integral of 
the dipole scattering amplitude Nh over the impact pa- 
rameter b, 



°dipv 



(r,a;) = 2 / dhN h {v,x;h). 



(2) 



dY J 2tt 

+ N(r 2 , Y) - N{r, Y) - N(r u Y)N{r 2 ,Y)\ , (3) 

where we define r = x — y, r\ = x — z, r 2 = y — z. Here, 
x (y) is the position of the q (q) in transverse space with 
respect to the centre of the target and z is the correspond- 
ing one for the emitted gluon. The BFKL kernel is given 

by 



K{r,r 1 ,r 2 ) = a s 



a,N c 



(4) 



In this framework, the QCD evolution is included in the 
dipole forward scattering amplitude. The simplest equa- 
tion describing this evolution and taking into account sat- 
uration effects is the Balitsky-Kovchegov (BK) equation 
[2]. We will first review the properties of the solutions of 
the BK equations in configuration space and then discuss 
whether these properties are observed in presently avail- 
able data. 



* Joint contribution to the proceedings of the Hard Probes 
2004 Conference based on the oral presentations by J.G. Mi- 
lhano and C.A. Salgado 



The BK equation © has a rather simple probabilis- 
tic interpretation [3]: when evolved in rapidity, the par- 
ent dipole with ends located at x and y emits a gluon, 
which, in the large- N c limit, corresponds to two dipoles 
with ends (x,z) and (z,y), respectively. The probability 
of such emission is given by the BFKL kernel ijljl. and 
weighted by the scattering probability of the new dipoles 
minus the scattering probability of the parent dipole (as 
the variation with rapidity of the latter is computed). 
The non-linear term, subtracted in order to avoid double 
counting, prevents, in contrast to BFKL, the amplitude 
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from growing boundlessly with rapidity. The BK equation 
ensures unitarity locally in transverse configuration space, 
l^( r 7^)l < 1- This is guaranteed since for N(r,Y) = 1, 
the derivative with respect to Y in © cannot be positive. 

The BK equation © was derived at leading order in 
a s In (s/so) for a fixed coupling constant a s . It is expected 
that next-to-leading-log corrections will play a significant 
role. An important part of these corrections will arise, as 
in BFKL, from the running of the coupling. The scale of 
the running coupling can only be determined in earnest 
once the full next-to-leading-log calculation is available. 
It is not clear a priori which of the three distance scales 
in the kernel l@J — an 'external' one, the size of the par- 
ent dipole; r, and two 'internal' ones, the sizes of the two 
newly created dipoles, r% and — should drive the run- 
ning of the coupling. In order to access the sensitivity 
of the results to an heuristically introduced running cou- 
pling, different modifications of the kernel (J3J were con- 
sidered in [4]. Essentially, the different physical cases are 
accounted for by three modifications of (0J: Kl, where 
the parent dipole scale r is used to evaluate the running 
of the coupling; K2, where the sizes r\ and r-i of the cre- 
ated dipoles drive the running of the coupling; and K?> 
which further modifies K2 by exponentially weighting the 
gluon emission vertex, thus imposing short range inter- 
actions, in order to check for possible sensitivity of the 
results to the Coulomb tails of the kernel. The coupling 
has been allowed to run in the standard one-loop form 
with three flavours. It has been frozen in the infrared to a 
value a s (k = 0) = ao, see [4] for details. 

In our numerical implementation of BK evolution (for 
a detailed discussion see [4]) we consider three different 
initial conditions evolved from some fixed value of xq (in 
practice, xq ~ 0.01). 

First, we consider an initial condition with the same 
r-dependence at fixed x$ as the Golec-Biernat-Wiisthoff 
model (GBW) [5], albeit with the x-dependence given by 
the BK evolution 1 



N' 



GBW 



(r) = 1 — exp 



(5) 



The second initial condition takes the form given by the 
McLerran-Venugopalan model [6] (MV): 



N MV (r) = 1 - exp 



r 2 Q' s 



In 



1 



r 2A2 
/l QCD 



(6) 



For transverse momenta k ~ 1/r > O(lGeV), the sen- 
sitivity to the infrared cut off e is negligible. The am- 
plitudes N GBW and N MV are similar for momenta of 
order Q' s but differ strongly in their high-fc behaviour. 
The corresponding unintegrated gluon distribution 4>(k) = 



1 Here and in the other initial conditions ©, below, we 
denote by Q' 3 what is usually called the saturation scale. Our 
definition of the saturation scale Q 3 is somewhat different, see 
Equation @ below, but the relation between both scales is 
straightforward e.g. in GBW, Q? = -4 In (1 - k) Ql- 



I 4^ eir ' k N ( r ) decays exponentially for N GBW but has 
a power-law tail ~ 1/k 2 for N M . Finally, we consider 



A^ s (r) = l-exp[-(rQ' s ) c ] 



(7) 



The interest in this ansatz is that the small-r behaviour 
N AS cx r c corresponds to an anomalous dimension 1 — 7 = 
1 — c/2 of the unintegrated gluon distribution at large 
transverse momentum. This anomalous dimension can be 
chosen to differ significantly from that of the initial condi- 
tions N GBW and N MV . Our choices c = 1.17 and c = 0.84 
are somewhat arbitrary. They can be motivated a poste- 
riori by the observation that the anomalous dimension of 
the evolved BK solution for both fixed and running cou- 
pling lies between the anomalous dimension of the initial 
conditions N AS and N GBW (or N MV ). Thus, the choice 
of N AS is very convenient to establish generic properties 
of the solution of the BK equation. 



3 Numerical results 
3.1 Evolution 

Figure^shows the evolution of the dipole scattering prob- 
ability for GBW initial condition with fixed and running 
coupling. The evolution is much faster for fixed coupling 
than for running coupling. Further, the differences arising 
from the specific implementation of running coupling ef- 
fects in the modified BFKL kernels Kl, K2 and K2> are 
rather small and well understood qualitatively [4] . Impor- 
tantly, effects of imposing short range interactions, K3, 
are very small (unless the range of the interaction is un- 
physically small). 



N{r) . Running, K1 

1 _ Fixed, KO 




. K1 


Y=0,6,12,18 


K2 




. --■ K3 




N GBW (Y=0) / 









Fig. 1. Solutions of the BK equation for GBW initial condition 
(dotted line) for rapidities Y = 6, 12 and 18 with Qo = 0.4. 
Left plot: Evolution with fixed (K0, solid lines) and running 
coupling (Kl, dashed lines). Right plot: evolution with running 
coupling for kernel modifications Kl (solid lines), K2 (dashed 
lines) and K3 (dashed-dotted lines). 
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3.2 Geometrical scaling 

In the limit Y — > oo, the solutions of the BK evolution are 
no longer functions of the variables r and Y separately, 
but instead depend on a single scaling variable 



rQ s (Y). 



(8) 



The saturation momentum Q S (Y) determines the trans- 
verse momentum below which the unintegrated gluon dis- 
tribution is saturated. It can be characterized by the po- 
sition of the falloff in N(r), e.g. via the definition 



N{r = l/Q S (Y),Y) = K , 



(9) 



where n is a constant which is smaller than, but of order, 
one. Different choices such as K = 1/2 (as in the results 
given below) and k = 1/e lead to negligible differences in 
the determination of Q S (Y). 

The accuracy of scaling at small r can be checked by 
comparison with the analytical scaling forms, for fixed and 
running coupling respectively, proposed in [7] : 



f)( T ) = aT 2 -> (lnr 2 +6) 



(10) 



| Fixed coupling u =0.2 

□ N 3 ™ 

g *N* S (C=1.17) 

A N AS (c=0.84) 
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Fig. 2. The rapidity dependence of the parameter 7, charac- 
terizing the anomalous dimension 1 — 7, as determined by a 
fit of < lot to the BK solutions for different initial conditions: 
GBW (squares), MV (circles), and AS with c = 1.17 (stars) 
and c = 0.84 (triangles). Left plot: results for fixed coupling 
with arj = 0.2. Right plot: results for running coupling with 
a<j = 0.4 and two versions of the kernel Kl (empty symbols) 
and K2 (filled symbols). 



r 27 h nT 2 + ± 



(11) 



Here, 1 — 7 is usually called the anomalous dimension 
which governs the leading large-A: behaviour of the un- 
integrated gluon distribution. 

We determine 7 from a fit of our numerical results to 
the functions (|10[) and (|llf> in the Y-independent region 
1(T 5 < t < 1CT 1 , i.e. for 10 5 Q s > 1/r > 10 Q s , with a, 
7 and S as free parameters. The results given below were 
found to be insensitive to a variation of the lower limit of 
this fitting range. 

For the case of fixed coupling constant, we find that the 
function Z 1 ) provides a very good fit to the evolved solu- 
tions. In Figure[21 we show the fit values of the parameter 
7, obtained for fixed coupling constant from the evolution 
of different initial conditions N GBW , N MV , and N AS for 
different values of c. At initial rapidity, these distributions 
have widely different anomalous dimensions but evolution 
drives them to a common value, 7 ~ 0.65, which lies close 
to the theoretically conjectured one [7, 8] of 0.628. For 
a small fixed coupling constant a?o = 0.2, this asymptotic 
behaviour is reached at Y ~ 70, while for a larger coupling 
constant <5o = 0.4 the approach to this asymptotic value 
takes half the length of evolution (results not shown) . For 
fixed coupling solutions, f 2 ^ does not provide a good fit 
to our numerical results. 

We have repeated this comparison for all running cou- 
pling prescriptions. We found that both f 1 ^ and f 2 ' pro- 
vide good fits and yield very similar values of 7. The re- 
sults for K3 are numerically indistinguishable from those 
for K2 and will not be shown in what follows. Also, the 
value of 7 was found to be independent of the coupling 
constant o.q at r — » 00. In Figure [21 we show the values 
of 7 extracted from a fit to f x K Irrespective of the ini- 
tial condition, they approach a common asymptotic value 



7 ~ 0.85. While our numerical findings for N AS with 
c = 0.84 are not inconsistent with the approach to this 
asymptotic value, no firm conclusions can be drawn. This 
initial condition just starts too far away from the asymp- 
totic scaling solution to reach it within the numerically 
accessible rapidity range. In this case, the monotonic in- 
crease of 7 with rapidity at large Y is smaller than the 
increase for N AS with c = 1.17 at comparable values of 
7, indicating that the rapidity evolution of the anomalous 
dimension depends in general not only on the small-r be- 
haviour, but on the full shape of the scattering probability. 

The value 7 ~ 0.85 is considerably larger than the 
one found in fixed coupling evolution. This is in agree- 
ment with previous numerical results [9] but in contrast to 
theoretical expectations [7, 8, 10] which predict the same 
value of 7 for the fixed and running coupling cases. As 
an additional check, we have performed running coupling 
evolution from an initial condition given by the solution 
at large rapidity of fixed coupling evolution (for which 
7 ~ 0.65). We find that even with this initial condition, 
running coupling evolution leads to a value of 7 ~ 0.85. 

It has been argued [7,8] that expressions ljl0|l and (|1 l|l 
are only valid for values of r inside the scaling window, 
t sw ~ Aqcd/Qs{Y) < t < 1 with Yq the initial rapidity, 
and that the dipole scattering probability returns to the 
doublc-leading-log (DLL) expression 

A^ DiL (r) = o(Y) r 2 [- In (r 2 yl 2 )]- 3 / 4 x 



x exp b(Y) y/- In (r 2 A 2 ) , (12) 

with a(Y) oc Y 1 / 4 and 6(F) ex y/Y, for values r < r sw . We 
have checked that this form provides a good fit (fit and 
numerical solution differ by less than ±10%) to the fixed 
coupling solution of BK for r < r sw = A/Q S (Y), A ~ 0.2 
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where a s = <5o — constant, A — dao and Qq = Q 2 (Y = 0) 
(i.e., the evolution starts at Y — 0). 

For running coupling, the momentum scale is expected 
to be ~ Qs(Y), suggesting 2 the substitution a s — ► a s (Q s (Y)) 
in Equation ipHfl) . This leads to [8] 



Q 2 S {Y) = A 2 exp 



AWY + X 



(16) 



where (A') 2 = 24N c d/f3 and X = {A')- 2 \n{Q 2 / A 2 ). 

The Y"-dependence of Q 2 for several initial conditions 
and for different choices of ao, calculated for all the kernels 
considered in this work, is shown in Figure 0] 



Fig. 3. Plot on the left: solutions of the BK equation (solid 
lines) with GBW initial condition and fixed coupling a 3 = 0.2 
compared to fits (dashed lines) to the DLL expression (12) i 
for rapidities Y = 0, 4, 8, 12, 16 and 20 (curves from right 
to left). Plots on the right: values of the coefficients a(Y) and 
b(Y) (circles) in the DLL expression versus Y , compared to fits 
(curves) to the functional form suggested by DLL. 



GeV, see Figure |3] Our comparison is limited to rapidi- 
ties Y < 20, since the scaling window starts to extend 
over the entire numerically accessible r-space for Y > 20. 
Up to Y = 20, the coefficients a(Y) and 6(F) follow the 
expected DLL V-behaviour, see Figure [3J However, the 
scaling ansatz f 1 ^ provides an equally good fit to the BK 
solutions for t < r sw . This is the reason why in previous 
numerical studies [11] no upper bound for a scaling win- 
dow was found. When the solutions of BK are fitted to 
Z 1 ) within the scaling window, the values of 7 at Y = 
for both initial conditions are < 20% smaller than those 
found when the fit is done within a fixed r-window. But 
for larger Y the values of 7 extracted from fits within ei- 
ther the scaling window or some fixed r-window approach 
each other and quickly coincide. 



3.3 Rapidity dependence of the saturation scale 

In the scaling region and for large rapidity (where Q S (Y) ^> 
Aqcd), the rapidity dependence of the saturation scale is 
determined [8] by 



d In [Q 2 (Y)/A 2 
dY 

where the numerical value of 
d 2 rd 2 Ti 1 



= dt 



(13) 



2vr 2 



'/= / " '.'V -T^[A r (ri) + ^(T 2 )-jV(r)-iV(T 1 )iV(T 2 )] 

can only be found once the scaling solution N(t) is known. 

For a fixed coupling constant, the saturation scale grows 
exponentially with rapidity 
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Fig. 4. The rapidity dependence of the saturation momentum 
Q 2 3 for fixed a s — 0.4 (thick solid), fixed a s — 0.2 (thin solid), 
and running coupling with ao = 0.4 for kernels Kl (dashed), 
K2 (dashed-dotted) and K3 (dotted lines). For each group, 
lines from top to bottom in the rightmost side correspond to 
initial conditions AS with c = 1.17, MV and GBW. 



The rise of Q 2 is much faster for fixed than running 
coupling, in accordance with l|15|) and l|16(l and as already 
observed in [9, 10, 14-19] 

For fixed coupling constant, Q 2 exhibits with good 
accuracy an exponential behaviour for high-enough val- 
ues of Y. The value of the slope extracted from a fit 
to the function l|15fl is A ~ 1.83 for ao = 0.4. As ex- 
pected, for ao = 0.2 this value is reduced by a factor two, 
A ~ 0.91. For the constant O, we find d ~ 4.57, in 
agreement with previous numerical studies at very high 
rapidities [11] but slightly smaller than the theoretical ex- 
pectation d = 4.88 [7,8]. 

For the case of a running coupling constant, an expo- 
nential fit can be done only for a very limited K-region. 
For example, for Y ~ 10 we find a logarithmic slope ~ 0.28 
for GBW or MV initial conditions with Q ~ 1 GeV, in 
agreement with the results of [10] but smaller than the val- 
ues found in [18] (see also [17,19]). The exponential func- 
tion (|15|) is unable to fit the full V-range. In contrast, the 
weaker rapidity dependence of l|16fl does provide a good 
fit in the full V-range. The fit to yields A' ~ 3.2, 



Ql{Y)=Ql exp [AY], 



(15) 



This approximation is also supported by numerical results 
[12-14] which show that in momentum space the typical trans- 
verse momentum of the gluons is ~ Q s - 
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while the theoretical expectation [7, 8] is slightly larger, 
A' = 3.6. 



where h Q 2 (Y = 0) is the initial saturation momentum for 
the nucleus, and (A 1 ) 2 is defined below Equation (|16fl . 



3.4 Nuclear size dependence of the saturation scale 

The nuclear size enters the initial condition. Here, we ex- 
amine the effect of BK evolution on the initial nuclear size 
dependence of the saturation scale. Figure [5] summarises 
our results. 



- A 1/3 


Y=0,10,40,72 ******** 


... A Mc ,c=1.17 










Fig. 5. Upper plot: QIa/QIp versus A 1 ^ 3 for initial conditions 
GBW (Q 2 sA {Y = 0) oc A 1/3 , solid) and AS with c = 1.17 
(Q 2 sa(Y = 0) oc A 2/3c , dashed lines); thick lines are the results 
for Y = in the running coupling case and for all rapidities in 
fixed coupling; for running coupling, different rapidities Y = 
10, 40 and 72 (thin lines) are shown from top to bottom for 
each initial condition. Lower plot: Q^aHA 1 Q%) versus Y 
for GBW with A 1/3 = 3, 6, 10, 20 and 50 with the same line 
convention as the upper plot (the results for fixed coupling 
have been obtained for Y < 36 and extrapolated as a constant 
equal to 1). In all plots ao = 0.4 and in the running coupling 
case the kernel Kl has been used. 



In the fixed coupling case, the initial ^-dependence of 
the saturation scale is preserved, irrespectively of whether 
this dependence is oc A 1 / 3 as for the GBW or MV ini- 
tial conditions (which produces numerical results for the 
A-dependence which are very close to those obtained for 
GBW), or whether it differs from oc A 1 / 3 due to an anoma- 
lous dimension as the one included in the AS initial condi- 
tion. This is in agreement with the theoretical expectation 
that, given the dilatation invariance of the BK equation 
J3J, any nuclear dependence present in the initial condi- 
tion can be scaled out. 

The evolution with running coupling is seen to reduce 
the A-dependence with increasing rapidity. If fitted in a 
wide rapidity range, the dependence of In [Q1a(Y)/Q 2 (Y)\ 
on Y is ~ Y~ aA . However, for large values of A and Y, 
the decrease with increasing Y is oc 1/y/Y and thus well 
described by [20] 



In 2 



In 



hQ 2 s (Y=0) 



QIa(Y) _ 



(17) 



4 Geometric scaling and experimental data 

4.1 Geometric scaling in lepton-hadron collisions 

Following [21] and motivated by the previous sections, let 
us assume that both the energy and the nuclear size (or 
centrality) dependence of the scattering amplitude N(r, x; b) 
can be encoded in the saturation scale Q s ,h(x, b) for any 
h (proton or nucleus). Then, the cross section can be 
written as 



(18) 



In this case, since L (Q 2 , r, z)\ 2 is proportional to Q 2 
times a function of r 2 Q 2 , the cross section is only a func- 
tion of t 2 = Q 2 /Q 2 h {x). This geometric scaling was found 
to describe all lepton-proton data with x < 0.01 [22]. In 
order to compare protons and different nuclei [21] one 
needs to make some assumptions about the impact pa- 
rameter dependence. Here, we have assumed that all the 
b dependence can be scaled by the nuclear radius of the 
hadronic target 3 h,b = b/^irR 2 , with R A = (1.12A 1 / 3 - 
0.86^4 -1 / 3 ) fm. Then, the condition for geometric scaling 
in lepton-nucleus data is 



al2(x,Q 2 ) = nRi j dvJdz\^l L (Q 2 1 v 1 z)Y 
x 2 [ dbN h (rQ s , h {x,b)). 



t1 *a 



(ta) 



* P {TA) 



nR 2 A 



ttR 2 p 



(19) 



For the A-dependence, we make the ansatz that the sat- 
uration scale in the nucleus grows with the ratio of the 
transverse parton densities to some power 1/6, which we 
take as a free parameter 



) 2 



( a^r 2 ; 



A-kR? v 



(20) 



Here, ttR 2 is the second free parameter to be fixed by the 
data. 

For the proton case we take the Golec-Biernat and 
Wusthoff (GBW) parametrization [5], Q 2 p = (x/xq)~ x in 
GeV 2 , x = 3.04 • 10~ 4 and A = 0.288. This parametriza- 
tion shows geometric scaling as can be seen in Fig. El In 
order to proceed to the nuclear case, we need a functional 
form of the scaling curve. The data [23] are seen to be well 
parametrized by [21] 

a^P(x,Q 2 ) =<P(t 2 ) = a Q [ lE + r (0,0 +ln£] , (21) 



This is exact for trivial impact parameter dependences as 
a step-function or a gaussian. We have checked that it gives a 
rather good approximation for a realistic, Wood-Saxon, profile. 
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where 7e is the Euler constant, r (0, £) the incomplete r 
function and £ = a/r 2b , with a = 1.868 and b = 0.746. 
The normalization is fixed by cto = 40.56 mb. 

To determine Q 2 A , we use Eqs. (tT§j) and lf2"UI) and com- 
pare the functional shape (|21[1 to the available experimen- 
tal data for -f*A collisions with x < 0.0175 [24-26], using 
£ = a/rf. We obtain S = 0.79±0.02 and 7ri? 2 = 1.55±0.02 
fm 2 for a x 2 /dof = 0.95 (see Fig. Elfor the comparison). 
Notice that these parameters translate into a growth of 
the saturation scale faster than A 1 / 3 for large nuclei. If 
we impose Q 2 A ~ ^l 1 ^ 3 in the fit, a much worse value of 
X 2 /dof = 2.35 is obtained. We conclude that the small- a; 
experimental data on 7* A collisions favours an increase of 
Q 2 A faster than A 1 / 3 . The numerical coincidence b ~ S is 

consistent with the absence of shadowing in nuclear par- 

2 

s,A- 



ton distributions at Q 2 3> Q 2 



3. 




J3 
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Fig. 6. Geometric scaling for 7*p (upper panel, data from [23]), 
"/*A (middle panel, data from [24,25]) and the ratio of data for 
"/"A over the prediction from 111' II (lower panel). As an addi- 
tional check, the lower plot also shows data for -y*A normalized 
with respect to 7*C [26], and divided by the corresponding pre- 
diction from Eq. H2H . 



4.2 Geometric scaling and multiplicities in AA 
collisions 

In [21] a simple formula for multiplicities in symmetric 
colliding systems at central rapidities has been proposed: 



1 dN 



AA 



N 



part 



dfj 



part 



(22) 



An easy way to derive this formula is to take, as the start- 
ing point, the factorized expression [27] for gluon produc- 
tion 



AD 



dN, 
dydptdh 



oc -§ dk (t) A (y,k 2 ,b)(f>B (y,{k- p t ) 2 ,b) , 
Pt J 

(23) 

where cf) h (y, k, b) = J dr exp{ir • k} N h (r, x; b)/(27rr 2 ) [28] 
If one assumes geometric scaling for the parton distri- 
butions, c/>A(y, k 2 ,b)= 4>(k 2 /Q 2 A (y, b)), we obtain 



dN. 



AA 



dy 



QIa 



Pf 



dkdh 



s,A 



(k-Pt) 



^drdb 0(r 2 )0((r 



(24) 



where the integrand is a constant. This proportionality 
between the total multiplicities and the saturation scale is 
shared by other models of hadroproduction [27-30]. It is 
important to notice, however, that for the case of geomet- 
ric scaling, Eq. i|22[l is more general than the factorized 
form l|23fl . Indeed, any functional shape of the integrand 
will lead to the same result providing geometric scaling 
holds. In order to recover Eq. (|22p. the energy dependence 
of the saturation scale in l|24() is given by the GBW param- 
eter A = 0.288; for its centrality dependence we use the 
known proportionality in symmetric A+A collisions be- 
tween the number -/V part of participant nucleons and the 
nuclear size A and Q 2 A cx A 1 / 35 , with S = 0.79 ± 0.02. 
In this way, the energy and centrality dependences are 
determined by parameters fitted to 7* — p and 7* — A re- 
spectively. In all these models, the hadron yield is assumed 
to be proportional to the yield of produced partons. The 
remaining normalization constant is fixed to No = 0.47 in 
order to reproduce PHOBOS data, see Fig. It is inter- 
esting that even the p+p data ( [31], as quoted in [32]) 
at y/s = 19.6 and 200 GeV are accounted for. Eq. 
implies that the energy and the centrality dependence of 
the multiplicity factorize, in agreement with the results by 
PHOBOS [32]. 

A formula similar to Eq. I|22|l has been extensively em- 
ployed to study multiplicities in Au+ Au collisions in [33] . 
These authors assume Q 2 A oc A 1 / 3 and an additional en- 
hancement factor ln(y / s A A f p^ t ) argued to come from scal- 
ing violations of the coupling constant in 123fl . Notice that 
for the accessible range of A, A i/9 — A 1 / 3 ln(A 1/3 ) - this 
being the reason why both approaches provide a fair de- 
scription of the data at RIIIC. 



4.3 Geometric scaling and dAu data 

The forward rapidity region of the RHIC experiments has 
become the main testing ground for saturation ideas. In 
this section, we check to which extent geometric scal- 
ing can account for the observed suppression on particle 
yields. The situation in this case is, however, more model- 
dependent than in the previous two cases. The reason is 
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Fig. 7. Energy and centrality dependence of the multiplic- 
ity of charged particles in Au+Au collisions 1221 compared to 
PHOBOS data [32] . Also shown in the lower panel are the p+p 
data [31] and results for y/a = 62.5 and 5500 GeV/A. 



the existence of two different saturation scales (one for 
the deuteron and another one for the gold nucleus) which 
precludes the trivial changes of variables done before. In 
order to proceed, we observe that if one writes Eq. in 
momentum space k, the main contribution to the cross sec- 
tion comes from k ~ Q/2. Approximating the dipole wave 
function (in momentum space) by 8{k — Q/2), the scal- 
ing curve (|21(l is (except for a normalization constant) the 
unintegrated gluon distribution, 4>a (k = Q/2) ~ <P(ta), 

with ta = fc 2 /4Q Sj A ■ For the case of particle produc- 

- 2 

tion in dAu collisions, the gluon saturation scale Q s ,a~ = 
N c Ql a /Cf needs to be used. Now, if the parton distri- 
bution in the deuteron falls off sufficiently quickly, <f>d ~ 
l/fc t ™, n > 1, we obtain from Eq. l(2"3]) 



dN° 



N coni dr,d2p t _ iVcoll 2 <MPt/Qs,c 1 ) _ N co n 2 $(T cl ) 
" N coni <j) A {pt/Q s ,c 2 ) ~ N coUi $(t C2 ) 



dA rdAu 
N coU2 dr]d 2 p t 



(25) 



for the centrality classes Ci, c-i. This simple formula relates 
the suppression measured in lepton-nucleus collisions to 
that in d+ Au. For the comparison in Fig. [S] to data [34] 
on the normalized yields of central and semi-central over 
peripheral d+Au collisions, we use the number of nucleon- 
nuclcon collisions AT co j] in different centrality bins [34] with 
AT colli = 13.6±0.3, 7.9+0.4 and iV coll2 = 3.3±0.4. Only the 
two most forward rapidities r\ — 2.2 and 3.2 are compared. 
This simplistic analysis indicates that the suppression of 
particle yields in forward rapidity measured in d+Au colli- 
sions is in agreement, through geometric scaling, with the 
nuclear shadowing measured in lepton-nucleus collisions. 
The connection between the small x- and A-dependence 
of parton distribution functions, and the suppression of 
normalized yields in d+Au collisions [34] at forward ra- 
pidity has been discussed in several recent works [11,35]. 
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Fig. 8. Normalized ratios of central and semi-central to periph- 
eral d+Au collisions measured by BRAHMS [34] compared to 
results from Eq. 12511 . The bands represent the uncertainty in 
the determination of JV C oii [34]. 



Eq. H25[) contributes to this discussion by illustrating to 
what extent the suppression of high-p t particles in d+Au 
at RHIC can be accounted for by the shadowing in 7* A 
collisions, see Eq. 1)21(1. 



5 Comparison of data and theory and 
conclusions 

The two main properties of the BK evolution are the ex- 
istence of a scale which indicates the presence of a sat- 
urated region and the existence of a scaling solution at 
asymptotic energies. The analytic form of the asymptotic 
solution is only known for some ranges of the dipole size. 
Here we have presented a numerical analysis of the BK 
equations without impact parameter dependence. In order 
to study the possible effects of next-to-leading-log correc- 
tions to the BK equation we have included the running of 
the coupling constant, using several prescriptions. Our re- 
sults confirm the existence of a scaling solution also for the 
running coupling case. Most of our numerical results agree 
with the analytical estimations. Namely, the evolution in 
rapidity of the saturation scale has been found to be ex- 
ponential for the case of fixed coupling and slower (the ex- 
ponential of a square root) for the running coupling case 
and with parameters in agreement with analytical esti- 
mates; the nuclear size dependence of the saturation scale 
has been found to be preserved in the fixed coupling case 
and to disappear asymptotically in the running coupling 
case; the small-r behaviour of the scaling solutions repro- 
duces the anomalous dimension found analytically for the 
case of a fixed coupling, but it is found to be different (in 
disagreement with analytical estimations) in the running 
coupling case. 

If the dynamics of the BK equations is relevant for 
the present experimentally accessible regimes, the above 
properties should be realized in the data. The observa- 
tion [22] that all the available data on lepton-proton scat- 
tering with x < 0.01 can be described by a single variable 
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t 2 = Q 2 /Ql p is reminiscent of the existence of a scaling 
solution. The energy dependence of the saturation scale 
extracted from a fit to experimental data [5] is too small 
compared with the results of BK with fixed coupling. It is, 
however, in approximate agreement with the running cou- 
pling case if one is restricted to a small range of energies. 
In the nuclear case, we have found [21] that a similar geo- 
metric scaling of the nuclear data is possible by encoding 
the nuclear size in the saturation scale. The A-dependence 
is, however, different from the usually assumed A 1 / 3 and 
turns out to be stronger, ~ A 4 / 9 . If this result is not a 
numerical accident, this behaviour could only have a dy- 
namical origin (the geometrical behaviour of Q 2 A is A 1 / 3 ). 
This is not in disagreement with the solutions we have ob- 
tained from the BK equations, but the relation is obscure: 
the yl-dependence of the saturation scale is preserved (i.e. 
it is the same as in the initial conditions) for the fixed cou- 
pling case, but disappears for the running coupling case. 
The nuclear size dependence is, however, expected to be 
modified when an appropriate treatment of the impact 
parameter is taken into account in the BK equation [36]. 

We have also found that a simple explanation for the 
centrality dependence of the multiplicities measured in 
symmetric systems at central rapidities is possible with 
the nuclear size dependence of the saturation scale ob- 
tained from lepton-nucleus data. In the presence of a ge- 
ometric scaling, these multiplicities are proportional to 
Qs a- This translates into a proportionality of the mul- 
tiplicities with the number of participants for the geomet- 
rical estimation Q 2 A ~ A 1 / 3 . The ^-dependence Q 2 A ~ 

^4/9 

gives, however, a nice description of the experimental 
data. As a final check, we have found a striking similarity 
of the suppression of nuclear yields in forward d+Au col- 
lisions measured by the BRAHMS collaboration at RHIC 
and the nuclear shadowing measured in lepton-nucleus ex- 
periments. All these observations are in agreement with 
the interpretation of the data in terms of saturation physics, 
but the existence of a numerical accident cannot be ex- 
cluded. A quantitative description of the experimental data 
with explicit use of QCD evolution equations will be needed 
in order to establish the relevance of the saturation effects. 

J. L. A. and J. G. M. acknowledge financial support by the 
Ministerio de Educacion y Ciencia of Spain (grant no. AP2001- 
3333) and the Fundacao para a Ciencia e a Tecnologia of Por- 
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